Pascal Super Library/Pascal Super Library (CW International)(1997).bin/MATH/PI64KPAS/PI_64K.PAS

Pascal/Delphi Source File | 1992-11-23 | 14.6 KB | 522 lines

{
      A Program to compute Pi upto 65531 didgets (about 64k didgets)

                  "PI_64K.PAS"  by Mitchell Ross  [70324,241]  11-23-92
Based on          "PIMUSIC.PAS" by Cino Hilliard  [72756,672]  10/10/87
that was based on "pi2.pas"     by Chuck White    [75006,3677] May 4, 1986
that was based on "pi1.pas"     by Cino Hilliard  [72756,672]  April 6, 1986
that was based on "pi.pas"      by Cino Hilliard  [72756,672]  ???

-----------------------------------------------------------------------------
                 Original Docs by Cino Hilliard Follow....
-----------------------------------------------------------------------------
this program computes the digits of pi using the arctangent formula

                (  term 1  )   (  term 2  )
    (1)  pi/4 := 4 arctan 1/5 - arctan 1/239

                                   in conjunction with the gregory series
                     n
    (2)  arctan x := sum [ (-1)^m*(2m + 1)^-1*x^(2m+1) ]  (approximately)
                    m:=0

substituting into (1) and (2) the first few values of n, rearranging,
simplyifing, and nesting  we have,

pi:= 3.2 + 1/25[-3.2/3 + 1/25[3.2/5 + 1/25[-3.2/7 + ...].].]       ( term 1 )
    -1/239[4 +1/239^2[-4/3 +1/239^2[4/5 +1/239^2[-4/7 +...].].]   ( term 2 )

using the long division algorithm and some tricks i discovered, this
( nested ) infinite series can be used to calculate pi to a large
number of decimal places in a reasonable amount of time.  a time
function is included to show how slow things get when n is large.

improvements can be made by increasing the number of digits each array
element holds and by changing the data type from integer to real for
selected variables. of course the added number of digits these changes
produce will cost much much more time. ah indeed, 'tis no free lunch!
however, since term 1 and term 2 are computed seperately and since the
arrays are step by step updated, the  program does lend itself to
parallel or non von ( what a coincidence - see below )  processing.

for example, let computer 1 perform term 1 and computer 2 perform  term 2.
moreover, let several computers share in the calculation of each
of the individual terms. however, to keep each computer equally busy,
a logarithmic type of adjustment must be made to decide on the number
of terms to be assigned to each computer. since the higher power terms
require much longer to compute, the computers assigned to higher powers
must be given  fewer terms to do.

a little history
----------------
in august, 1949, professor john von neumann used formulas  (1) and (2)
to calculate pi to 2035 decimal places on the  eniac  computer.
the effort was made to determine if the digits conformed to some type of
pattern or if they were randomly distributed.
the calculation was completed over the labor day weekend with the
combined efforts of four eniac staff members working in eight-hour
shifts to ensure continuous operation of the eniac. the calculation
(including card handling time) took approximately 70 hours.

the conclusion was as suspected - the digits appeared to be random!

some years ago i requested information on pi from the encyclopedia
britannica research service. i received a report giving the above
historical account plus a listing of the 2035 digits.

a little comment
----------------
using my t2k, pi1.pas computes pi to 2035 places in  15 min  31.97 sec.

gee whiz!  what earthly use can be made of this?  i plan to make designs
and color patterns using the digit values as the basis. maybe something
like fractals would be a good start. this is why computation and speed
are important as i do not wish (nor trust) to key in entries from a
published list.

                        cino hilliard
                        [72756,672]
-----------------------------------------------------------------------------
                       End of original Docs...
-----------------------------------------------------------------------------

My version in not in assembly, so that data types can easily by changed. The
old program would only calc to a little over 2 thousand, so I decided to see
just how many I could do. This program is the best I could do without going
to a random disk file (S..L..O..W..!!!!) or to EMS (alas, beyond me). Because
of the large compute times, I added a save/restore feature. Running this
program at night, durring work hours, and on weekends on my home 386-25 took
about 2 weeks to do the full 65531 didgets. I shudder to think what disk-based
arrays would have done to the compute time. Well, here it is! If anyone can
find a way to do more didgets, and still keep the speed up, please let me know.

             - Mitchell Ross
               3670 Road 5-1
               Delta, OH  43515      CI$: [70324,241]
}

{$I-}       { no file IO      (speed...) }
{$R-}       { no range checks (speed...) }
{$S-}       { no stack checks (and more speed!) }
{$N-} {$E-} { no co-processor, it's all integer math here }

Program CalcPi;

Uses CRT;

Const ArraySize = (64*1024)-2;         { as big as we can go... }
      OutFileName = 'pi.!!!';         { final output file }

      SaveFileName1 = 'PI_WORK.$$1';      { save files }
      SaveFileName2 = 'PI_WORK.$$2';

Type PiType = Array[0..ArraySize] of ShortInt;
     PiPtrType = ^PiType;

var HowFar : Integer;
    NumberToCalc : LongInt;                 { num of didgets }
    j, q, v, r : LongInt;
    Steps, Steps1, Steps2 : LongInt;
    ArrayMax : LongInt;

    p,t,a : PiPtrType;               { here's the 3 BigNum arrays }

{---------------------------}
{ Handy to check for existing files. Don't leave home without it. }
Function Exist( FileName : String ): Boolean;

Var WorkFile : Text;

Begin
 {$I-}
 Assign( WorkFile, FileName );
 Reset( WorkFile );
 If IOResult <> 0 then Exist := FALSE
  ELSE begin
   Exist := TRUE;
  Close( WorkFile );
 end;
end;

{---------------------------}
{ saves _everything_ to disk, so we can pick up again later. }
Procedure SaveSoFar;

Var SaveFile1 : Text;
    SaveFile2 : File of PiType;

begin
  Writeln;
  Write('Saving...');

  Assign( SaveFile1, SaveFileName1 );
  Rewrite( SaveFile1 );
  Writeln( SaveFile1, ArrayMax );
  Writeln( SaveFile1, HowFar );
  Writeln( SaveFile1, NumberToCalc );
  Writeln( SaveFile1, Steps );
  Writeln( SaveFile1, Steps1 );
  Writeln( SaveFile1, Steps2 );
  Writeln( SaveFile1, j );
  Writeln( SaveFile1, q );
  Writeln( SaveFile1, v );
  Writeln( SaveFile1, r );
  Close( SaveFile1 );

  Assign( SaveFile2, SaveFileName2 );
  Rewrite( SaveFile2 );
  Write( SaveFile2, a^ );
  Write( SaveFile2, p^ );
  Write( SaveFile2, t^ );
  Close( SaveFile2 );
  Writeln('Done!');
end;

{---------------------------}
{ and load it back again.... }
Procedure RestoreSoFar;

Var SaveFile1 : Text;
    SaveFile2 : File of PiType;

begin
  Writeln;
  Write('Restoring...');

  Assign( SaveFile1, SaveFileName1 );
  Reset( SaveFile1 );
  Readln( SaveFile1, ArrayMax );
  Readln( SaveFile1, HowFar );
  Readln( SaveFile1, NumberToCalc );
  Readln( SaveFile1, Steps );
  Readln( SaveFile1, Steps1 );
  Readln( SaveFile1, Steps2 );
  Readln( SaveFile1, j );
  Readln( SaveFile1, q );
  Readln( SaveFile1, v );
  Readln( SaveFile1, r );
  Close( SaveFile1 );

  Assign( SaveFile2, SaveFileName2 );
  Reset( SaveFile2 );
  Read( SaveFile2, a^ );
  Read( SaveFile2, p^ );
  Read( SaveFile2, t^ );
  Close( SaveFile2 );
  Writeln('Done!');
end;

{---------------------------}
Procedure CheckForQuit;

Var Key : Char;

begin
  Key := '▒';         { my favorite char, but just make it <> #27 }

  While KeyPressed do   Key := ReadKey;          { clear buffer }

  If Key = #27 then begin            { if escape pressed, save }
    SaveSoFar;
    HALT;
  end;
end;

{---------------------------}         { real guts of the program follow... }
{---------------------------}
procedure div_32_by_Steps_to_a;  { divide 3.2 by Steps and store in array a }

begin
 q := (3 div Steps);
 r := (3 mod Steps);
 a^[0]:= q;
 v := r * 10 + 2;
 q := (v div Steps);
 r := (v mod Steps);
 a^[1] := q;

 for j:= 2 to ArrayMax do begin
   v := r * 10;
   q := (v div Steps);
   r := (v mod Steps);
   a^[j] := q;
 end;
end;

{---------------------------}
procedure diva(d:integer); { divide a by specified integer  }

begin
  r := 0;
  for j := 0 to ArrayMax do begin
    v := r * 10 + a^[j];
    q := (v div d);
    r := (v mod d);
    a^[j] := q;
  end;
end;

{---------------------------}
procedure div_4_by_Steps_into_a;   { divide 4 by Steps and store in a }

begin
  q := (4 div Steps);
  r := (4 mod Steps);
  a^[0] := q;
  for j := 2 to ArrayMax do begin
    v := r * 10;
    q := (v div Steps);
    r := (v mod Steps);
    a^[j] := q;
  end;
end;

{---------------------------}
procedure div_32_by_Steps_into_p;  { divide 3.2 by Steps and store in p }

begin
  q := (3 div Steps);
  r := (3 mod Steps);
  p^[0] := q;
  v := r * 10 + 2;
  q := (v div Steps);
  r := (v mod Steps);
  p^[1] := q;
  for j := 2 to ArrayMax do begin
    v := r * 10;
    q := (v div Steps);
    r := (v mod Steps);
    p^[j] := q;
  end;
end;

{---------------------------}
procedure div_4_by_Steps_into_p;   { divide 4 by Steps and store in p }

begin
  q := (4 div Steps);
  r := (4 mod Steps);
  p^[0] := q;
  for j := 1 to ArrayMax do begin
    v := r * 10;
    q := (v div Steps);
    r := (v mod Steps);
    p^[j] := q;
  end;
end;

{---------------------------}
procedure suba;         { subtract a from p  and store in a }

begin
  for j := ArrayMax downto 0 do begin
    a^[j] := p^[j] - a^[j];
  end;
end;

{---------------------------}
procedure sub32a;        { subtract a from 3.2 and store in t }

begin
  t^[0] := 3;
  t^[1] := 1;
  for j := ArrayMax downto 2 do begin
    t^[j] := 9-a^[j];
  end;
end;

{---------------------------}
procedure sub4a;        { subtract a from 4 and store in a }

begin
  a^[0] := 3;

  for j := ArrayMax downto 1 do begin
    a^[j] := 9-a^[j];
  end;
end;

{---------------------------}
procedure subt;         { subtract term2 from term 1 and store in a }

begin
  for j := ArrayMax downto 1 do begin
    a^[j] := t^[j] - a^[j];

    while a^[j] < 0 do begin
      a^[j] := a^[j] + 10;
      a^[j-1] := a^[j-1] + 1;
    end;

    while a^[j] > 9 do begin
      a^[j] := a^[j] - 10;
      a^[j-1] := a^[j-1] - 1;
    end;

  end;
end;

{---------------------------}
procedure compute_a1;

begin
  CheckForQuit;               { give user chance to save work }
  Steps := Steps1;           { compute term 1 to Steps1 div 2  places }
  div_32_by_Steps_to_a;
  diva(25);
  HowFar := 1;
  CheckForQuit;               { give user chance to save work }
end;

{---------------------------}
procedure compute_a2;

begin
  while Steps > 3 do begin
    CheckForQuit;               { give user chance to save work }
    GotoXY(1,WhereY-1);
    Writeln( 'Pass 1:  ', Steps,' of ',Steps1,'  ');   { show progress }
    Steps := Steps -2;
    div_32_by_Steps_into_p;
    suba;
    diva(25);
  end;
  HowFar := 2;
  CheckForQuit;               { give user chance to save work }
end;

{---------------------------}
procedure compute_a3;

begin
  sub32a;
  HowFar := 3;
  CheckForQuit;               { give user chance to save work }
end;

{---------------------------}
procedure compute_b1;

begin
  Writeln;
  Steps := Steps2;          { compute term 2 to Steps2 div 2 places }
  div_4_by_Steps_into_a;
  diva(239);
  diva(239);
  HowFar := 4;
  CheckForQuit;               { give user chance to save work }
end;

{---------------------------}
procedure compute_b2;

begin
  while Steps > 3 do begin
    CheckForQuit;               { give user chance to save work }
    GotoXY(1,WhereY-1);
    Writeln(  'Pass 2:  ', Steps,' of ',Steps2,'  ');  { show progress }
    Steps := Steps -2;
    div_4_by_Steps_into_p;
    suba;
    diva(239);
    diva(239);
  end;
  HowFar := 5;
  CheckForQuit;               { give user chance to save work }
end;

{---------------------------}
procedure compute;      { compute the nested series for Steps1 iterations }

begin
  If HowFar <= 0 then  compute_a1;          { if we did it before, don't }
  If HowFar <= 1 then  compute_a2;          { do it again.               }
  If HowFar <= 2 then  compute_a3;

  If HowFar <= 3 then  compute_b1;
  If HowFar <= 4 then  compute_b2;

  sub4a;
  diva(239);
  subt;                 { subtract term 2 from term 1 and store in a }
end;

{---------------------------}
procedure printpi;    { print the  formated digits of pi }

Var OutFile : Text;
    z : Integer;

begin
  Writeln('░▒▓█ DONE! Saving.... █▓▒░');
  Assign( OutFile, OutFileName );
  Rewrite( OutFile );

  write( OutFile, 'π = 3.');
  for j:=1 to NumberToCalc do begin
    write( OutFile, a^[j]);
    if j mod 5 = 0 then write( OutFile, ' ');
    if j mod 50= 0 then writeln( OutFile, ' ',j:4,' pl');
  end;
  writeln( OutFile );
  writeln;
  writeln('  Results saved to ',OutFileName );
  writeln;

  Close( OutFile );

  {$I-}
  Assign( OutFile, SaveFileName1 );          { done now, kill temp files }
  Erase( OutFile );
  Assign( OutFile, SaveFileName2 );
  Erase( OutFile );
  z := IOResult;          { throw away any errors }
end;

{---------------------------}
procedure start;               { prompt for input and initialize }

begin
   HowFar := 0;

   If Exist( SaveFileName1 ) then begin
     Writeln('------ Restoring from saved session -----');
     RestoreSoFar;
   end ELSE begin
     write('Number decimal places? (25..',ArraySize-3,') :');
     readln(NumberToCalc);
     Writeln;

     ArrayMax  := NumberToCalc+2;
     Steps1 := 4*(3*ArrayMax div 8)+3;   { number of terms- series 1 }
     Steps2 := 4*(Steps1 div 14)+3;      { number of terms- series 2 }
   end;
  Writeln;
  Writeln('Press ESC to save work so far. (Might take a moment after keypress)');
end;

{-----------------------}
{-----------------------}

var z : Integer;

begin
  ClrScr;
  Writeln('╔═════════════════════════════════════════════════╗');
  Writeln('║ Pi_64k - By Mitchell Ross                       ║');
  Writeln('║             3670 Road 5-1                       ║');
  Writeln('║             Delta, OH  43515   CI$: [70324,241] ║');
  Writeln('║                                                 ║');
  Writeln('║   Compute the value of PI, up to 64k didgets    ║');
  Writeln('╚═════════════════════════════════════════════════╝');
  Writeln;

  New( a );
  New( p );
  New( t );

  start;
  compute;
  printpi;
end.